% Jay Angel

% Progam similating a neuron through circuits

%%%%%%%%%%%%%%%%%%% Variables %%%%%%%%%%%%%%%%%%
%% ELS Comments
% Error 1 The  b vector  needed  negative signs for the chosen orien
%
%Error 2The time vector, at 1 seconsd was too short to see the assymptote
%
%Error 3 The y plotting range needed to allow for negative swings.

%% ELS General Comments
% Some intution about the expected asymptotic voltage and time constant
% should have been used to avoid some of those errors.
%
% But, good job, almost there!

% A = Reduced incidence matrix 
% b = Voltage source vector
% bNa = Sodium voltage source
% bK = Potassium voltage source
% G Conductance matrix
% gNa Sodium resistor
% gK Potassium resistor
% C Capacitance matrix 
% e potential differences
% x potential at node 1


%%%%%%%%%%%%%% Derivation of Essential Equations %%%%%%%%%%

% G = -C d\dt e = -C d\dt(b-Ax)
% dx\dt = (At*C*a)^-1*[At*G*A-At*b*A]
% dx = (At*C*a)^-1*[At*G*A-At*b*A]*dt

%%%%%%%%%%%%%%%%%% Function %%%%%%%%%%%%%%%%%%%%
clc
clear
close all

%%Incidence matrix
% A= -1 1;1 -1;-1 1
% A= [-1 1 1]^t

%Creates the Incidence Matrix
i=[];j=[];s=[]; 
i=[i 1]; j=[j 1]; s=[s -1];
i=[i 1]; j=[j 2]; s=[s 1];
i=[i 2]; j=[j 1]; s=[s 1];
i=[i 2]; j=[j 2]; s=[s -1];
i=[i 3]; j=[j 1]; s=[s -1];
i=[i 3]; j=[j 2]; s=[s 1];
A=sparse(i,j,s);


%% Reduced incidence matrix
%This removes the last column from A simulating the ground in the circuit
A=A(:,1); 


%% Conductivity matrix and Capacitance Matrix
%.5 micro Semens
gNa=.5*10^-6;
%5 micro Semens
gK=10*10^-6; 

G=[gNa 0 0;0 gK 0;0 0 0];
G=sparse(G);
  
%If the neuron is assumed to be a sphere of radius 20 um with a specific
% capacitence of 1uF/cm^2 the resulting capacitance is
%els 
%els C=(pi*(20*10^-6)^2)*(100)^2; % used capital C on left and right hand side of eqn
% els lost factor of 4 in surface area c=(pi*(20*10^-6)^2)*(100)^2;
%els you left out the specific capacitance
s_c = 1e-6; % els micro Farads/cm^2
c = 4*pi*(20*10^-6)^2*s_c*1e4;
% els C=[0 0 0;0 0 0;0 0 C]; % used capital C on both sides!
C=[0 0 0;0 0 0;0 0 c];
G=sparse(G);
C=sparse(C);

 
%% Voltage sources
% els both of these should be negative
% These batteries are edge traversed anti-physically
bNa=-.055; % 55 mV source
bK=-.075; % 75 mV source
b=[bNa;bK;0]; % b is a column vector
b=sparse(b);

s=struct('G',G,'C',C,'b',b);
%% Function that uses incidence matrix A and the parameters saved in struct
%els next line needs to over more time than 100 clicks!
 x=zeros(100); % note this is solution vector 100 time points and a vector of size(nodes) as a matrix!
 t=zeros(100,1);
 %x(1)=0; t(1)=0;
 delt=1e-7;
 for j=2:length(x)
 x(j)=x(j-1)+delt*A'*s(1).G*(s(1).b-A*x(j-1))/(A'*s(1).C*A);
 t(j)=t(j-1)+delt;
 end
 % next line y limits needs to be negative
 %figure,plot(t,x); title('RC\_Strang'); axis([0 10 0 .1]);ylabel('Volts');xlabel('Time-Secs');
 figure,plot(t,x); title('RC\_Strang'); axis([0 1e-5 -.1 .1]);ylabel('Volts');xlabel('Time-Secs');
 title('Correct asymptote at -69mv; Correct time constant');
 %%Publish to html 
 %needs to be copied into another script in order to work but this makes it
 %easier to send it to you
  %publish('RCNeuron.m');
  
 
 
 